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ABSTRACT 

We present a full hydrodynamical investigation of a system of self-contained, interact- 
ing gas 'cloudlets'. Calculations are performed using SPH, and cloudlets are allowed 
to collide, dissipate kinetic energy, merge and undergo local gravitational collapse. 
The numerical feasibility of such a technique is explored, and the resolution require- 
ments examined. We examine the effect of the collision time and the velocity field 
on the timescale for evolution of the system, and show that almost all star formation 
will be confined to dense regions much smaller than the entire cloud. We also inves- 
tigate the possibility, discussed by various authors, that such a system could lead to 
a power-law IMF by a process of repeated cloudlet coagulation, finding that in fact 
the inter-cloudlet collisions occur at too high a Mach number for merging to play an 
important part. 

Key words: hydrodynamics - methods: numerical - stars: formation - ISM: clouds. 



1 INTRODUCTION 

Star forming regions are often interpreted in a 'cloud-clump- 
core' picture, in which an overall system (the cloud) consists 
of an ensemble of discrete regions of gas (clumps), which 
move about with supersonic velocities through some confin- 
ing medium. The collisions and mergers of these clumps lead 
to regions of gravitational collapse (cores) and eventually to 
star formation. This idea is largely an attempt to charac- 
terise the observed structure of molecular cloud complexes 
(see for example the review by Williams, Blitz & Mc Kee 
I2OOOI) , but is likely to be applicab le in many situations. 
Variou s authors such as lKwanI il979h an d ISilk fc Takahashil 
il979l) have cited coagulation p rocesses as an origin f or the 
form of the IMF. For example, iMurrav fc LinI Jl99d) have 
argued that in protoglobular clusters, gas is fragmented into 
sub- Jeans mass clumps by a variety of hydrodynamical in- 
stabilities, and that eventual star formation ensues when 
enough clumps have merged to form Jeans unstable objects. 
A similar process may be involved in the assembly of Gi- 
ant Mole cular Clouds (GM Cs) in the spiral arms of disc 
galaxies dPringle et aLllioOlll . since the diffuse ISM that is 
gathered together in galac tic shocks is inferred to be highly 
clumped on small scales jKovama fc Inutsukall200ol) . It is 
the interactions of these clumps, or 'cloudlets', that controls 
the formation of stars, and understanding the properties of 
such systems is important in order to shed light on questions 
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such as the timescales for star formation and the origin of 
the IMF. 

Until recently, a full treatment of such a system was 
beyond the computing resources available. Previous authors 
have examined numerically the interaction of gas cloudlets, 
but the majority of these studies have been limi ted to the 
collision of two objects. Early work JChieze fc Laz arcfi IQS^; 
lHausr nan"l981 ^ was severely limited by the resolution avail- 
able. iLattanzio fc HenriksenI (4988') carried out a thorough 
examination of the p ossible outcomes of a two-body col- 
lision. L ater studies ([ Pongracic et alJ Il992!: 'Bhatta l et alJ 
[1998; Ki mura fc TosalllOT 6: Mari nho fc Lep inc 2000) ah ex- 
amine in detail the structure produced in one collision be- 
tween two large bodies of gas. 

None of these numerical investigations concern a sys- 
tem of_jnan2_^nteracting bodies. One attempt to do 
so jMonaghan fc Varnas, .1983) examined a group of 48 
cloudlets using SPH, but was limited to a very low res- 
olution. Instead, authors have overcome the computing 
limitations by using adapted N-body methods, in which 
the results of cloudlet collisions are determined accord- 
ing to pres cribed rules whenever two cloudlets approach 
each other iScalo fc Pumphrevlll982l: iMurrav fc Linill99d 
iMurrav. Ravmondao^^^^^anskiT^ggl) . A full hvdrodv- 
namical treatment of such a system using Smoothed Par- 
ticle Hydrodynamics is only now becoming computationally 
feasible, and it is to this problem that we devote this paper. 

Our aim is not to present a realistic model of any as- 
tronomical object, but rather to investigate the numerical 



2 D. Gittins, C. Clarke and M. Bate 



modelling of such ensembles. We will examine the feasibil- 
ity of such simulations, with particular regard to the nu- 
merical resolution required, and discuss the limitations of 
this approach. We will also explore the general properties 
of such systems, with particular regard to the energy dis- 
sipation timescale, the distribution of protostars and the 
role of coagulation in the evolution of the system. In this 
way we w ill be able to r eassess earlier investigations such 
as that of lMurrav fc Lir] lIlQQSl . who found that successive 
clump-clump collisions built up a spectrum of resulting 'stel- 
lar masses' that was broadly compatible with a power law 
IMF. 

In addition to the abstract questions described above, 
this paper also lays the ground work for a series of papers in- 
vestigating the formation of GMCs from a shocked, clumpy 
ISM in the spiral arms of disc galaxies. A key feature of these 
latter simulations involves the behaviour of colliding clumps 
in the complex velocity field generated in spiral shocks, and 
thus it is of considerable importance for these more realistic 
simulations that the generic properties of coalescing clump 
ensembles, and the numerical requirements for such calcula- 
tions, are well understood. 

The structure of the paper is as follows: in section|21we 
set out the numerical method and the results of a suite of 
two-body clump-clump encounters. The aim here is to de- 
fine regions of parameter space associated with particular 
outcomes (merger, collapse, or clump destruction) and to 
define the number of particles per clump that is required 
in order to achieve numerical convergence of expected out- 
comes. Section 13 contains a detailed description of the sys- 
tem investigated and outlines the key parameters. In sec- 
tion 01 we describe generic properties of the cloud evolution 
and focus on the following aspects: the relationship between 
the energy dissipation timescale and the nominal two-body 
collision timescale ( thus re-evaluating the claims made by 
IScalo fc PumphrevI lll982h in this regard) , the possibility of 
such a system leading to a distributed cluster of protostars 
by collisionally-induced local gravitational collapse, and the 
process of building up a realistic IMF/clump mass spectrum 
by this process. Section |S] summarises our chief conclusions. 



The code makes use of two common modifications. The 
first concerns the fate of a region of self-gravitating mate- 
rial that begins to collapse. Unchecked, the material will fall 
inward towards infinite density, and as it does so the CPU 
time required to advance the gas rises without limit. This 
problem can be o vercome by the use of 'sink particles' (see 
iB^eeTai] mil) for more details). In this technique, when 
an SPH particle passes a critical density (pcrit = 10® times 
the initial mean density), a sink particle is formed by replac- 
ing the SPH gas particles contained within race = O.lrd oudlct 
of the densest gas particle in a pressure-supported fragment 
by a point mass with the same mass and momentum. Any 
gas that later falls within this radius is accreted by the point 
mass if it is bound and its specific angular momentum is less 
than that required to form a circular orbit at radius race from 
the sink particle. Sink particles interact with the gas only 
via gravity and accretion. Sink particle mergers are not al- 
lowed. Each sink particle therefore represents an area of the 
gas which has collapsed under its own gravity. The further 
fate of such a region is clearly not resolved by the calcula- 
tion. These sink particles will be referred to as 'protostars' 
in the remainder of the paper. 

The second modification is the inclusion of a confining 
external pressure. The cloudlets in our simulations are as- 
sumed to be moving through a thin confining medium, con- 
sisting of a hot gas that exerts sufficient pressure to keep the 
undisturbed cloudlets confined. The cloudlets do not other- 
wise interact with the confining medium, and exchange no 
momentum with it. The pressure of this external medium 
can be easily included in SPH by introducing a 'constant 
pressure' technique. When the pressure force between two 
particles is calculated, the pressure of each particle p is re- 
placed by a modified pressure p' = p — Pcxt, where poxt is 
the constant external pressure, fixed at the start of the cal- 
culation. 

The calculations presented here used between 10"" and 
lO'^ SPH particles. The largest of these require significant 
computing power to run in a practical length of time, and 
the UKAFF supercomputing facility in Leicester was used 
for this purpose. 



2 NUMERICAL TECHNIQUE 
2.1 Hydrodynamical code 

The calculation presented here was performed using a three- 
dimensional, smoothed particle hydrodynamics (SPH) code. 
The S PH code is based on a versio n originally developed by 
Benz jBendll990l : lBenz et alJll990t) . The smoothing lengths 
of particles are variable in time and space, subject to the con- 
straint that the number of neighbours for each particle must 
remain approximately constant at A^ncigh = 50. The SPH 
equations are integrated using a second-order Runge-Kutta- 
Fehlb erg integrator with individua l time steps for each par- 
ticle (iBate. Bonnell fc Pricelll995^ . Gravitational forces be- 
tween particles and a particle's nearest neighbours are cal- 
culated using a bi nary tree. We use the standard form of 
artifi cial viscosity jMonaghan fc Gingoldl Il983l : iMonaghaiil 
^93) with strength p arameters Qy = 1 and /3v = 2. Further 
details can be found in lBate et al.Nl995^ . The code has been 
parallelised by M. Bate using OpenMP. 



2.2 Resolution requirements 

The resolution requirements must be carefully determined to 
ensure that the calculation is always resolved, whilst keeping 
the total number of SPH particles to a minimum. This allows 
the calculations to be completed in a reasonable length of 
time while following the evolution of the system for as long 

as p ossible. 

iBate fc Burkerii il997l) determined that the basic re- 
quirement for an SPH calculation of this type is to ensure 
that the number of SPH particles per Jean's mass remains 
above 2Ancigh particles. In the simulations presented here, 
the cloudlets will collide with each other, and in the col- 
liding region an isothermal shock layer will be formed. In 
such shocks, the density is increased by a factor equal to 
the square of the Mach number M . The Jean's mass in such 
layers (and thus the number of particles per Jean's mass) 
will consequently fall by a factor M. In this way, given the 
expected Mach number of collisions, a sufficient number of 
particles can be chosen to ensure that the resolution require- 
ment will be met in the isothermal shock layers. 
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In addition, two more tests are required to demostrate 
that the appropriate hydrodynamical processes are being 
correctly simulated - firstly, the kinetic energy dissipation 
in a collision between two cloudlets, and secondly, the con- 
ditions under which such a collision can lead to a collapse 
or a merger. These two issues are discussed in the following 
sections. 
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2.2.1 Kinetic energy dissipation in cloudlet collisions 

In a collision between two cloudlets, kinetic energy will be 
dissipated. Since this is the dominant process in the evolu- 
tion of the systems we investigate, it is important to check 
that this loss of energy is followed sufficiently accurately at 
the numerical resolution we have used in our calculations. 

We checked this requirement by performing a single 
two-body encounter with a range of nine resolutions, from 
50 up to 5 X 10^ SPH particles per cloudlet. The calcula- 
tion follows the collision of two spherical cloudlets of ra- 
dius 1 at Mach 4, at an impact parameter b/r = 1 (i.e. the 
separation perpendicular to the direction of travel is equal 
to the radius of one sphere). The calculation is performed 
in dimensionless units such that the gravitational constant 
G, the sound speed of the gas Cs and the mass of each 
cloudlet are all equal to 1. The cloudlets are Bonnor-Ebert 
spheres with their bou ndary at ^ ~ 2.08 (using the termi- 
nology of iBonnoilllQSd) . corresponding to a density contrast 
of Pccntro/pouter — 1-8. The external pressure is chosen to 
confine the spheres at the appropriate radius. Under these 
conditions each sphere contains approximately 0.25 Jean's 
masses. The cloudlets are given initial velocities of u = 2 
toward each other. At time t = the separation along the 
direction of travel is 4 radii, such that if the cloudlets did not 
interact, their centres would pass each other at time t = 1. 

The progress of the collision is shown in figure Ini- 
tially, the kinetic energy increases as the cloudlets attract 
each other. When they hit, a shock layer forms. The pres- 
sure gradient tends to push particles away from the direction 
of travel, and within the shock layer, much of the mate- 
rial from one cloudlet avoids being reduced to zero velocity 
by 'slipping past' material from the other. As the cloudlets 
separate, material is spread out between them. As this gas 
expands, its reduced pressure will create a further deceler- 
ation on the trailing edges of the initially non-overlapping 
regions of the cloudlets, dissipating further kinetic energy. 
This expanding region develops a complex structure at high 
numerical resolution, which is unresolved when fewer parti- 
cles are used. The surviving portions of each cloudlet escape 
with a reduced velocity. The total kinetic energy of the gas 
is displayed in figure |5| 

As is clear from figure |21 the greater the numerical res- 
olution, the less kinetic energy is dissipated in the collision. 
This is to be expected, since at lower resolution the smooth- 
ing lengths of the particles will be necessarily larger, and 
consequently, any dissipative region will be 'smoothed' over 
a larger volume, and hence a larger proportion of the mass. 
In figure 13 the kinetic energy of the system at time t = 4 is 
plotted against the logarithm of the number of particles. 

It is somewhat difficult to assess the 'realism' of this 
result. A theoretical value for the proportion of kinetic en- 
ergy dissipated is not readily obtainable. Various previous 
authors have examined similar collisions, but few have con- 
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Figure 1. Progress of the two-cloudlet collision used to test the 
kinetic energy dissipation in collisions. Two cloudlets of isother- 
mal gas are collided at Mach 4. The cloudlets are Bonnor-Ebert 
spheres with g ~ 2.08. The initial separation along the direction 
of travel was 4 radii, and perpendicular to the direction of travel 
was 1 radius. The upper four panels show the calculation using 
50,000 particles per cloudlet; the lower four, using 100 particles 
per cloudlet. 



sidered the energy dissipated; where such a value has been 
quoted, the effect of altering the numerical resolution has 
not been examined. It would appear from figure|2lthat there 
may be a turnover around A'^part = lO'', implying that the 
calculation would converge on some value at ever increas- 
ing resolution. It is clear that using only 100 particles per 
cloudlet will overestimate the amount of energy dissipated 
in each collision, but the magnitude of this overestimation is 
not clear. We conclude that this is a fundamental limitation 
of attempting to model molecular clouds using this tech- 
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Figure 2. Kinetic energy vs. time in tlie collision shown in figure 
m The units are dimensionless (the gravitational constant, sound 
speed, initial radii and mass of each sphere are 1). In the absence 
of any interaction, the centres of the cloudlets would pass each 
other at time t = 1. Nino different resolutions are shown, with ap- 
proximately the following numbers of SPH particles per cloudlet: 
(top to bottom) 500,000, 50,000, 15,000, 7,500, 3000, 1000, 500, 
100 and 50. 
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Figure 3. Kinetic energy at time t = 4 in the two-cloudlet col- 
lision shown in figure plotted against log(A'^) where N is the 
number of SPH particles per cloudlet. 



nique, and that one should therefore avoid relying on the 
detailed numerical value of the kinetic energy in such sys- 
tems, exploring instead the overall evolution and structure 
of the cloud. 



2.2.2 Merging and collapse of colliding cloudlets 

The general problem of the collision of two clouds of gas 
has been examined by many authors, as discussed in section 
The range of conditions that can be simulated is very 
large, an d the range of outcomes is e qually great (see for 
example iLattanzio fc Henriksenlll98?j) . Any analysis must 
be restricted to a subset of the parameter space. We intend 
here to examine the effect of varying numerical resolution on 
the outcome of a simple two-body collision, of the type that 
dominates the calculations presented later in this paper — 
namely, the collision of two identical Bonnor-Ebert spheres. 

There are many possible outcomes from the collision 
of two such cloudlets. Work on simulating such encounters 
shows that these outcomes can be divided into three cat- 
egories: a) a collapsing object, in which the central region 
becomes gravitationally unstable, b) a merger, in which the 
end result is one large, stable cloudlet containing all the mass 
of the two progenitors, or c) a 'splash', in which instance the 
overlapping portion of the gas tends to be ripped from the 
cloudlets and splashed into the region between them, and 
two smaller regions of gas survive the encounter and escape 
with a reduced velocity in either direction. These three pro- 
cesses are central to the evolution of the systems described 
later in this paper. We carried out a range of calculations 
to determine the conditions under which each of these out- 
comes can be expected. Collisions were simulated over three 
initial values of the Mach number M , the Jeans fraction per 
cloudlet /j , and the impact parameter b. The result of each 
calculation was categorised as a collapse, merger or splash, 
as described above. Finally, to ensure that the minimum 
resolution requirements were met, the calculations were re- 
peated using 100, 300 and 1000 particles per cloudlet. 

The results of the calculations are shown in table 
They indicate that at low Mach numbers, the result is al- 
ways a merger or a collapse. The deciding parameter be- 
tween these outcomes is, as expected, the total number of 
Jeans masses, with the calculations at /j < 0.5 always pro- 
ducing a merger and those at /j > 0.5 always collapsing. In 
the borderline cases with /,j = 0.5, the more 'glancing' col- 
lisions (i.e. those at higher impact parameter) are less likely 
to collapse. At high mach number, the result is always a 
splash unless the collision is head-on. 

The results at the three different numerical resolutions 
were unchanged except in two 'borderline' cases, indicated in 
table0 In each of these calculations the result was a merger 
when 100 particles were used, and a collapse at higher reso- 
lution. This result is attributable to the number of particles 
involved in the overlapping region of the collision - at low 
resolution, there are simply too few particles in the merged 
region to create a sink particle. 

This investigation can be extended, for example to col- 
lisions between cloudlets of unequal size and/or mass, ro- 
tating cloudlets and so forth. However, our purpose here is 
not to undertake a full analysis of two-body collisions, but 
to investigate the numerical implications of the kind of col- 
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Table 1. Results of two-cloudlet collisions for a range of Jean's 
fraction (/j), Mach number {A4) and impact parameter (fe). Each 
collision ends in collapse (*), merger (o) or splash (-). The radial 
parameter ^ of the relevant Bonnor-Ebert density profile is also 
shown. Calculations marked f ended in collapse when > 100 par- 
ticles per cloudlet were used, and in merger when 300 or 1000 par- 
ticles were used. The remaining calculations produced the same 
outcome at all three numerical resolutions. 
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lisions that dominate the calculations presented later in this 
paper. This point is discussed further in sectional 



3 DESCRIPTION OF THE SYSTEM 
3.1 General description 

The system under investigation consists of many identical, 
spherical clouds of isothermal gas, which will hereafter be 
referred to as 'cloudlets'. These are pressure-confined by a 
tenuous surrounding medium, with which they do not oth- 
erwise interact. Each cloudlet initially has a stable Bonnor- 
Ebert density profile, and the external pressure is chosen to 
match this density profile at the chosen cloudlet mass and 
radius. The cloudlets are initially uniformly distributed in 
a spherical volume, and given individual bulk velocities in 
random directions. In most of the calculations the veloci- 
ties are of equal magnitudes. The speeds of the cloudlets 
are then scaled to fix the virial ratio Q = \Ek/Eg\ to the 
desired value (where Ek and Eg are the overall kinetic and 
gravitational energies of the cloud). 

These initial conditions are not intended to be a real- 
istic representation of a molecular cloud. Our calculations 
do not include magnetic fields, stellar feedback, heating and 
cooling of the gas, or any kind of 'support' mechanisms. 
This allows us to isolate the evolution of the system under 
the processes of gravity and hydrodynamics. The object is 
to investigate a system of gas cloudlets set up in such a way 
as to encourage the development of a cluster through the 
processes of collisions and mergers. 

In most of the calculations presented here, the cloudlets' 
velocities are scaled such that Q = 0.5. This is done in order 
to encourage the most interactions between cloudlets during 
the lifetime of the cloud. In super-virial clouds, the cloudlets 
will tend to fly apart and the mean free path will increase; in 
sub-virial clouds, cloudlets will simply fall to the centre. In 
section [4. II we demonstrate this point by examining clouds 
with increased and decreased values of Q. 

An isothermal equation of state has been used. This is 



a good approximation to the equation of state of molecu- 
lar gas at the densities likely to be relevant in applications 
of the models presented here. In molecular clouds, heating 
due to collisions and collapse radiates o n a shor t timescale, 
such that the gas is roughly isothermal llLarsoni ri969'l. This 
approximation holds in collap sing molecular clou d cores up 
to densities of ~ 10^^ g cm~^ llMasunaga fc Inut suka 2000)). 
Hence, an isothermal equation of state is an appropriate sim- 
plification for the simulation of gas systems such as those 
presented in this paper. 

3.2 Parameters 

The properties of this system are controlled by a few key 
parameters: 

(1) R, the overall radius of the cloud 

(2) A'^, the number of cloudlets 

(3) Cs, the sound speed of the gas 

(4) /j, the fraction of a Jean's mass initally contained in 
each cloudlet 

(5) r, the ratio of the collision timescale to the crossing 
time of the cloud (= icoii/icross), where taoss = 2i?/u and 
the cloudlets initially have a velocity v 

(6) Q, the ratio of kinetic to gravitational energy. 

Since each cloudlet is a pressure-bounded Bonnor-Ebert 
sphere, it will obey 

2 Gm 

Cs = « 

r 

where m and r are the mass and radius of the cloudlet re- 
spectively, and 

2 ^ _2 

Hence, choosing fj also fixes m/r for each fragment. Fixing 
r is equivalent to fixing X/R, the ratio of the mean free path 
to the radius of the cloud. Since 

3iVr2 

choosing X/R also fixes r/R for each fragment. 

Consequently, once the overall scaling is established by 
setting R and Cs, it is sufficient to select values of N, /j 
and r to determine the mass and radius of each cloudlet. 
The cloudlets must then be assigned bulk velocities, usu- 
ally with a random direction and all of the same magnitude 
(other prescriptions for the initial velocity distribution are 
described below). The velocities are then scaled to give the 
desired ratio Q of kinetic and gravitational energy. In most 
models, the cloud is virialised so that Q — 0.5. Finally, an 
important property of the system is the initial Mach num- 
ber of the cloudlets A4 . Since M indicates the typical Mach 
number expected in collisions, it is important in determining 
the minimum resolution, as discussed in section 

In all the calculations presented here, N = 1000. Table 
|5| lists all the models described in the text and the values of 
r, /j and Q chosen, as well as the velocity distribution used. 

We have simulated clouds in which all cloudlets have 
identical masses and radii. This is not an attempt to model 
(for example) a real molecular cloud, in which a range of 
'clumps' of different sizes would be seen. Rather, the inten- 
tion is to see whether a system with a range of cloudlet sizes 
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Figure 4. Mass and radius of stable Bonnor-Ebert spheres over 
a range of Jeans fraction /j. The temperature is fixed such that 
Cs = 1, and the external pressure fixed such that a sphere of mass 
m = 1 will have a radius r = 1. The units are dimensionless such 
that G = 1. Note that there is a maximum sphere radius, which 
occurs when Cg = Gm/2r. 



can arise entirely through the processes of two-cloudlet inter- 
ac tions, without being; s pecified at the outset (as suggested 
bv lMurrav fc LinI (|l99fih . for example). 

Furthermore, the possibility of setting up a range of 
cloudlet sizes in one cloud is limited by the restriction that 
the cloudlets are stable Bonnor-Ebert spheres. This is indi- 
cated in figure^ which shows the mass and radius of stable 
Bonnor-Ebert spheres at a range of /j, subject to a fixed 
temperature and external pressure. The stable solutions all 
have roughly the same radius, except at very low /j (at 
which values, a very large number of mergers would be re- 
quired to form a gravitationally unstable object), and so it 
would not be feasible to set up a cloud in which the cloudlets 
had a significant range of radii. Models with high values of 
/j are described in section \A/2\ 



4 CALCULATION RESULTS 

The 'baseline' calculation (model 1) used A'' = 1000 
cloudlets, /j = 0.25 and t = 1. Under these conditions, 
the cloud initially contains 250 Jean's masses and the colli- 
sion timescale tcoii is equal to the crossing time tcross. The 
virialisation of the cloud results in each cloudlet being given 
a bulk velocity of ~ 2.7. Consequently, random collisions 
between cloudlets can be expected to take place at Mach 
numbers up to around = 5. Each cloudlet contained 100 
SPH particles. Several random realisations of the same sys- 
tem were run to ensure that the result was independent of 
the specific initial conditions. 

The general evolution of the system is shown in figure|S| 
As the cloudlets collide, their relative kinetic energy is dis- 
sipated, and material is 'splashed' into the region between 
them as they separate. This material falls toward the centre 



Table 2. Definitions of the parameters used in the various models 



Velocities 
i 
i 
i 
d 
i 



/j = mcloudlet/A'-fj and Q = 
I^k/^gI are defined in the text. The letters in the Velocities col- 
umn refer to the initial velocity dispersion and are as follows: (i) 
isotropic and uniform, (d) isotropic with a distribution of speeds, 
(e) elliptically distributed and (c) spatially correlated. 
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Figure 5. Column-density maps of the cloud used in model 1, at 
times t = 0, t = 0.5icross, t = tcross and t = 1.25tcross. Protostars 
arc indicated by star symbols. 



of the gravitational well. After a time tcoii, most of the gas 
has formed a large dense region in the centre of the cloud. 
The gas in this region continues to collapse, and goes on 
to form a complex structure in which locally gravitation- 
ally unstable regions form. However, these processes are not 
resolved in our calculations, and consequently do not bear 
close examination. Such complex turbulent regions must cur- 
rently be studied separately using very high resolution (see 
for example Bate, BormcU & Bromm 2002). 

The calculation consists of two phases. To begin with, 
kinetic energy is dissipated approximately linearly with time 
as the cloudlets collide with one another. After each colli- 
sion, splashed material, which tends to lose specific angular 
momentum in the collision, begins to fall to the centre. This 
continues until roughly t = tcoii, at which point the average 
number of collisions per cloudlet is expected to be about 1. 
By this time, a significant fraction of the gas has fallen to 
the centre of the cloud. Cloudlets passing through this cen- 
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Figure 6. The decay of total energy vs. time in model 1 (t = 1, 
solid), model 2 (r = 0.2, short dashed) and model 3 (r = 5, long 
dashed) . 

tral region are disrupted and add to the gas already present. 
The system enters a second phase, in which the kinetic en- 
ergy is rapidly dissipated, and the dense gas at the centre 
begins to produce local gravitational collapses as described 
above. 

Within this dense central region, complex interactions 
of cloudlets (involving more than two cloudlets at a time) 
and disrupted material will occur. However, an examination 
of such processes is not the intent of this paper, and (as 
mentioned above) these regions are not resolved in our cal- 
culations. We will not, therefore, discuss here the further fate 
of material that has fallen into a region of runaway collapse, 
leaving such work to more detailed numerical investigations. 

Given this fact, the evolution of the system is dominated 
by two-body collisions. Interactions between more than two 
cloudlets (outside the central region) are rare. 

4.1 Dissipation timescale results 

It is of particular interest to examine the timescale on which 
the evolution of the system takes place. Figure El shows the 
change in the total (e.g. kinetic plus potential) energy of the 
system (in which the collision time tcou — Across) with time. 
It is clear that in this case, the kinetic energy is dissipated 
on the collision timescale, tcoii. Also shown in figure|^are the 
results from two further calculations with r = 0.2 and r = 5 
(models 2 and 3 respectively), in order to clearly indicate 
the dependence of the energy dissipation on icon- 

When tcoii is reduced by a factor of five, the initial rate 
of dissipation of kinetic energy is approximately five times 
greater as expected. The runaway dissipation begins at a 
time t ~ O.Stcross- This broadly agrees with the simple ex- 
pectation that cloudlets will begin to fall to the centre after 
a time of roughly tcoii (= 0.2tcross in this case), and take 
around half a crossing time to do so. Similarly, when tcoii is 



increased by a factor of five, the initial dissipation rate is re- 
duced by roughly the same factor. In this case, the runaway 
collapse phase is not reached by t = 2tc:ross. 

Scalo & Pumphrey (1982) performed calculations of a 
similar system of cloudlets, using a modified N-body code. 
The parameters of the system used in their calculation are 
essentially the same as that used in ours, except that they 
did not include gravity - the cloud was kept confined by us- 
ing a hard reflective outer boundary. They concluded that 
the dissipation time in such a system is at least an order of 
magnitude longer than the collision time, due to the propen- 
sity of collisions at low impact parameter. 

In our calculations, the system will not survive for more 
than roughly the collision time. There are two main reasons 
for this difference. The flrst is the 'splashing' of material 
in collisions in our calculations. The simulations of Scalo & 
Pumphrey used prescribed rules in which the outcome of any 
collision was either one, two or three discrete cloudlets, at 
the same density as the colliding clouds. Consequently, there 
were no spread-out regions of low density gas, which would 
otherwise have added to the total coUisional cross-section 
of the gas as the calculation proceeded. Secondly, since no 
gravity was included, the cloud did not relax toward a cen- 
tral density profile, and there was no preferred centre toward 
which gas would fall after a collision. Hence the formation 
of a central dense region, causing the runaway dissipation in 
the second phase of our calculations, was impossible. 

We also examined the effect of introducing a distribu- 
tion of initial cloudlet speeds, by repeating the calculation 
as follows (model 4). Cloudlets were given velocities in a 
random direction, and magnitudes chosen randomly from a 
probability distribution p{v), with 

p{v) oc exp(— u^). 

The velocities were then scaled as before to ensure overall 
viriality. All other parameters of the calculation were un- 
changed. In this case, the initial rate of energy dissipation 
was equal to that seen in model 1, but the cloud went into 
runaway collapse after around O.Stcross (compared to model 
1 which took twice as long to reach this phase). This differ- 
ence is attributable to the cloudlets initially at low speeds, 
which tend to simply fall to the centre of the cloud and ar- 
rive at the same time, creating a dense region and initiating 
the collapse. 

In two further runs, we examined the effect of altering 
the virial coefficient Q = \Ek/Eg\ (for all models except 
5 and 6, Q = 0.5). Model 5 is highly sub-virial, such that 
Q = 0.01, and is otherwise identical to model 1. In model 6 
the cloud is super- virial with Q — 1, and in order to satisfy 
resolution requirements, the Jeans fraction has been reduced 
to /j = 0.2. In these calculations the crossing time (as de- 
fined from the initial velocity of the cloudlets) is not very 
meaningful, so instead we will use the gravitational free-fall 
time, given by 



This is roughly half the crossing time tcross = 2R/v for viri- 
alised clouds, confirmed by setting «^ — GM/R in the above 
to recover 
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Figure 7. The decay of total energy vs. time in model 1 (Q = 0.5, 
solid), model 5 (Q = 0.01, long dashed) and model 6 (Q = 1, short 
dashed) . 
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The change in total energy in models 5 and 6 compared 
to model 1 is shown in figure |7| The results are largely what 
is expected. In model 5, the cloudlets simply fall to the cen- 
tre. Some energy dissipation occurs as they merge, but the 
dissipation rate is somewhat reduced compared to the viri- 
alised cloud. Several protostars are formed through these 
mergers, although by the time the first of these has formed, 
the maximum density has increased beyond the limit set 
by the resolution criterion set out in section 12.21 and con- 
sequently the gravitational collapse at all of these points is 
not guaranteed to be realistic. After about one free-fall time, 
the cloud goes into runaway collapse at the centre. In model 
6, however, the clouds undergo high Mach number collisions 
from the outset and the initial rate of energy dissipation is 
high. This rate gradually decreases as the cloud increases in 
size and the mean free path correspondingly grows. 



4.2 Protostar distribution results 

An important question to ask of the system under investiga- 
tion is whether it can lead to a cluster of protostars, and if 
so, what the spatial distribution of those protostars will be. 
These 'protostars' are represented in our numerical code by 
replacing collapsing regions of gas with 'sink particles'; see 
section r2. II for a description of the formation and accretion 
criteria used. 

In the baseline calculation, the resulting distribution of 
protostars was highly clustered in the centre of the cloud. All 
of these were formed by the fragmentation of the large mass 
at the centre - none were formed in two-cloudlet interactions. 
This result holds independently of the collision timescale of 
the cloud. 
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Figure 8. Column-density maps of model 7 (/j = 0.5) at time 
t = (left) and t = tcross (right). Protostars are indicated by star 
symbols. 



Figure 9. Column-density maps of model 8 (/j = 0.8) at time 
t = (left) and t = tcross (right). Protostars are indicated by star 
symbols. 



Since the Jean's fraction in each cloudlet, /j, is only 
0.25, four cloudlets need to be assembled together before 
one Jean's mass can be reached. If, however /j is increased, 
a smaller number of cloudlets will be required, and the for- 
mation of protostars in cloudlet interactions should become 
more likely. Hence, two further calculations were performed, 
with /j = 0.5 (model 7) and fj = 0.8 (model 8), in order 
to encourage as far as possible the formation of protostars 
throughout the volume of the cloud rather than in the centre 
only. 

The results of these calculations are shown in figures |S| 
and|5| Protostars are now formed in all parts of the cloud 
in two-cloudlet collisions, and as expected, more protostars 
are formed when /j — 0.8 than when /j — 0.5. However, 
at the end of the calculation, only a small proportion of the 
initial mass of the cloud (0.1 for model 7 and 0.15 for model 

8) is present in these protostars - the remainder once again 
forms a large mass at the centre of the cloud. Thus, although 
increasing /j significantly does result in the formation of 
some protostars throughout the cloud, the majority of the 
mass still falls into a small region in the centre. 

The central distribution of the mass in the final state of 
these calculations is a consequence of the isotropic velocity 
distribution. A more complicated initial velocity distribution 
might be expected to result in a corresponding change in the 
eventual outcome. In order to investigate this, further cal- 
culations were performed, in which the velocity distribution 
took the form of a prolate ellipsoid, such that the cloudlets' 
velocities were preferentially aligned along one axis (model 

9) . The results of these calculations are shown in figure 1101 

The effect of the velocity anisotropy is clear: rather than 
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Figure 10. Column-density maps of model 9 at time t = (left) 
and t = O.Stcross (right). Protostars arc indicated by star symbols. 



forming one large central mass, the gas has fallen into two 
dense regions at the foci of the velocity ellipsoid. Otherwise, 
the results are largely unchanged; the kinetic energy of the 
system is still dissipated in roughly 0.5fcoii, and local gravi- 
tational collapse is still confined to the two dense regions. 



4.3 Coagulation 

It is known from the previous results outlined in section 
I2.2.2l that cloudlets encountering each other at a sufficiently 
low Mach number (which is around 1) will merge to form 
a larger cloudlet. If a sufficient number of cloudlets merge 
together, at some point they will surpass the gravitational 
stability limit and start to collapse. If this process continues 
throughout the evolution of the cloud, the continual merg- 
ers will lead to a distribution of cloudlet mass, and many 
collapsed objects could be formed by this mechanism. 

This idea has be en discussed as an origin for the IMF in 
globular clusters bv lMurrav fc LinI iH-QQ&l . They examined 
somewhat similar systems to those discussed in this paper, 
using an N-body code in which cloudlets instantaneously 
merge when they touch. They found that this process re- 
sulted in a mass power spectrum of protostellar objects at 
the end of the calculation. It might, therefore, be expected 
that a similar process will occur in our simulations, leading 
to a power spectrum of final collapsed objects. 

In fact, very few mergers occur and the coagulation pro- 
cess is largely unimportant in our calculations. The main 
reason for this is simply that the great majority of collisions 
occur at too high a Mach number. The prescribed rule that 
any two cloudlets will merge on contact, used by Murray & 
Lin, is replaced in our work by the hydrodynamical result 
of the collision, which is (in general) a dissipative 'splash' 
rather than a merger. 

The typical Mach number of collisions in the system 
does, of course, depend on the initial parameters. If the ini- 
tial velocity distribution is isotropic and the velocities are 
scaled such that the cloud is virialised, then the relation be- 
tween the mach number and the remaining parameters A'', 
T ~ X/2R and a can be found. Starting from 

Ca 

and using v'^ ~ M/i?, Cg ~ am/r and M — Nm, we have 

--^^ 

Then, since 



N 



Using the parameters a — 1, N — 1000 and t — 1 the 
initial velocity of the clumps is around Mach 2.5, leading 
to collisions at up to Mach 5. If we wish to ensure that all 
collisions are subsonic, the initial velocities need to fall by 
a factor 5. For fixed and a, this implies r = 625. At 
this value, a given clump will be expected to have only one 
collision in 625 crossing times. Alternatively, if r is fixed 
then either (for fixed A*") a = 25, which corresponds to /j ~ 
2 X 10"^, or (for fixed a) N = 1.6. The latter case is clearly 
not interesting, and in the former, around half of all the 
cloudlets would need to merge before one Jeans mass could 
be reached. 

The conclusion is that the typical Mach number of colli- 
sions in the system, if the cloud is virialised and the velocity 
distribution is uniform, cannot be reduced to subsonic level 
without resorting to very long mean free paths, very stable 
initial cloudlets, or both. 

It is possible to encourage collisions at lower velocities 
by introducing a spat ially correlat ed initial velocity field. 
The Larson relations llLarsonl|l98 J) state that the velocity 
dispersion in star-forming clouds varies with the length scale 
I asl'^ , where a ~ 0.38. This correlation can be applied to an 
ensemble of cloudlets by repeatedly subdividing the volume 
of the calculation into eight cubes, then subdividing each 
cube again, and at each level assigning each cube random 
velocities on the appropriate scale. Murray & Lin examined 
systems both with and without this property. 

In the case of a cloud of 1000 cloudlets, subdividing the 
volume three times leaves only one or two cloudlets remain- 
ing in each subcube. Thus, the length scale at the smallest 
subdivision will be a factor 8 smaller than the overall cloud. 
The velocity dispersion on this scale will consquently be re- 
duced by a factor g"-^^ ~ 2.2 compared to the overall cloud. 
This can be expected to produce many more collisions at 
low Mach number, and hence promote the formation of pro- 
tostars due to agglomeration of several cloudlets. 

To this end, we repeated the calculation of model 1, 
with the introduction of a spatially correlated velocity field 
(model 10) . The cloud is repeatedly subdivided as described 
above, and the velocity for each subcube is chosen with a 
random direction and a magnitude picked from a gaussian of 
the appropriate width. The velocities are once again scaled 
to ensure overall viriality. The results are shown in figure 

EH 

In this calculation, mergers did occur, and several pro- 
tostars were formed outside the core of the cloud by the 
collapse of merged cloudlets. However, only a small num- 
ber of protostars were formed before the system once again 
formed a runaway collapse region at the centre. This oc- 
curred at f ~ 0.5fcross, as is expected, since the cloudlets in 
this model begin with a spread of initial speeds similiar to 
that used in model 4 (see section I^TTl . 

The mass distribution of cloudlets at this time is shown 
in figure 1121 The distribution of protostar masses is also 
shown, but note that the majority of these formed in the 
unresolved cloud centre with masses above m = 10 and 
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Figure 11. Column-density maps of model 10 at time t = 
(left) and t = O.Stcross (right). Protostars are indicated by star 
symbols. 
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linear) rate determined by the collision timescale, and a sec- 
ond collapse phase in which the majority of the gas falls 
to dense regions and becomes gravitationally unstable. This 
second phase will begin much sooner if a broad distribution 
of initial cloudlet speeds (rather than a 5-function) is used. 
The evolution timescale of the whole cloud is thus affected 
by the mean free path of the cloudlets, and the initial veloc- 
ity distribution also plays a significant part in determining 
the cloud lifetime, which can be significantly different from 
the dynamical timescale of the cloud. 

Protostar formation is mostly confined to the dense 
region that forms in the centre of the cloud, even if the 
cloudlets are initially close to gravitational instability. 

Finally, we find that the process of hierarchical assem- 
bly of cloudlets through mergers to form a mass spectrum 
of final objects is not significant. This is because the great 
majority of collisions occur at too high a Mach number and 
do not result in merging of the cloudlets. This will remain 
true so long as: the collision timescale is comparable to the 
crossing time of the cloud, and the cloudlets are not of such 
a low mass as to be very far from gravitational instability. 
Introducing a spatially correlated velocity field can promote 
mergers, but the effect is still small and few protostars will 
be formed by this mechanism. 
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Figure 12. Mass distribution of uncoUapsed objects (solid line) 
and protostars (dashed line) from model 10 at t = 0.5£cross. The 
vertical dotted lines indicate the initial cloudlet mass and initial 
Jeans mass. The slope of the Salpeter IMF is indicated by the 
long dashed line. 

are not shown. The peak at log(m) = represents undis- 
turbed cloudlets with m — 1, while the following two peaks 
correspond to m = 2 and m — 3, indicating the number 
of complete two-cloudlet and three- cloudlet mergers. Many 
cloudlets with masses below m = 1 are formed from material 
sheared off in collisions, and these are in fact more numer- 
ous than cloudlets formed at masses m > 2 from mergers, 
reinforcing the point that disruptive collisions play a more 
important role in the system than mergers. The slope of the 
Salpeter IMF (^abctcr 1955) is also shown, and it is clear 
that the distribution of merged cloudlets is too steep to be 
consistent with a stellar IMF. 



5 CONCLUSIONS 

The general evolution of a system of cloudlets as laid out 
in this paper tends to consist of two phases: an initial dis- 
sipative phase in which kinetic energy is lost at a (roughly 
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